Persiapan data¶

data sudah dilakukan cleaning sebelumnya

In [1]:
!pip install geopandas
Requirement already satisfied: geopandas in c:\users\user\anaconda3\lib\site-packages (0.13.2)
Requirement already satisfied: packaging in c:\users\user\anaconda3\lib\site-packages (from geopandas) (21.3)
Requirement already satisfied: pyproj>=3.0.1 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (3.6.0)
Requirement already satisfied: pandas>=1.1.0 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (1.4.2)
Requirement already satisfied: shapely>=1.7.1 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (2.0.1)
Requirement already satisfied: fiona>=1.8.19 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (1.9.4.post1)
Requirement already satisfied: six in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.16.0)
Requirement already satisfied: importlib-metadata in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (4.11.3)
Requirement already satisfied: cligj>=0.5 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: click~=8.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (8.0.4)
Requirement already satisfied: click-plugins>=1.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: attrs>=19.2.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (21.4.0)
Requirement already satisfied: certifi in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (2023.5.7)
Requirement already satisfied: colorama in c:\users\user\anaconda3\lib\site-packages (from click~=8.0->fiona>=1.8.19->geopandas) (0.4.4)
Requirement already satisfied: python-dateutil>=2.8.1 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2.8.2)
Requirement already satisfied: pytz>=2020.1 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2021.3)
Requirement already satisfied: numpy>=1.18.5 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (1.22.4)
Requirement already satisfied: zipp>=0.5 in c:\users\user\anaconda3\lib\site-packages (from importlib-metadata->fiona>=1.8.19->geopandas) (3.7.0)
Requirement already satisfied: pyparsing!=3.0.5,>=2.0.2 in c:\users\user\anaconda3\lib\site-packages (from packaging->geopandas) (3.0.4)
In [2]:
import geopandas as gpd
import warnings
warnings.filterwarnings(action='ignore')

loc_spring = gpd.read_file("Lokasi_mata_air.shp")

loc_spring.head(2)
Out[2]:
OID_ Name FolderPath BeginTime EndTime HasLabel LabelID kelas_pema luas_jangk Nama TDS__ppm_2 debitt2 lat long geometry
0 0 Bunut Document/Layer 2021-04-30 04:43:58.0 1 -1 0 NaN 0.0 Bunut 65 0.3 107.969 -6.83253 POINT Z (107.96867 -6.83253 0.00000)
1 0 Cikondang Document/Layer 2021-04-29 03:32:51.0 2 -1 0 NaN 3200.0 Cikondang 57 1.0 107.974 -6.84136 POINT Z (107.97449 -6.84136 0.00000)
In [3]:
loc_spring['EndTime'] = loc_spring['EndTime'].astype(int)
loc_spring.head(2)
Out[3]:
OID_ Name FolderPath BeginTime EndTime HasLabel LabelID kelas_pema luas_jangk Nama TDS__ppm_2 debitt2 lat long geometry
0 0 Bunut Document/Layer 2021-04-30 04:43:58.0 1 -1 0 NaN 0.0 Bunut 65 0.3 107.969 -6.83253 POINT Z (107.96867 -6.83253 0.00000)
1 0 Cikondang Document/Layer 2021-04-29 03:32:51.0 2 -1 0 NaN 3200.0 Cikondang 57 1.0 107.974 -6.84136 POINT Z (107.97449 -6.84136 0.00000)
In [4]:
import pandas as pd

# Baca file Excel (excel)
tambahan = pd.read_excel('mtairpython.xlsx')
tambahan.head(2)
Out[4]:
Kode Mata air PH Salinitas Suhu TDS Debit DHL warna Rasa Jarak patahan kelompok
0 1 Bunut 6.9 0 26.1 65.0 0.3 94 bening tawar 1.70 C
1 2 Cikondang 6.6 0 30.0 57.0 1.0 109 bening tawar agak sepat 1.04 D
In [5]:
# Gabung shapefile dan file Excel berdasarkan kolom ... dan ...
loc_spring2 = loc_spring.merge(tambahan, left_on='EndTime', right_on='Kode')

loc_spring2.columns
Out[5]:
Index(['OID_', 'Name', 'FolderPath', 'BeginTime', 'EndTime', 'HasLabel',
       'LabelID', 'kelas_pema', 'luas_jangk', 'Nama', 'TDS__ppm_2', 'debitt2',
       'lat', 'long', 'geometry', 'Kode', 'Mata air', 'PH', 'Salinitas',
       'Suhu', 'TDS', 'Debit', 'DHL', 'warna', 'Rasa', 'Jarak patahan',
       'kelompok'],
      dtype='object')
In [6]:
# Daftar kolom yang ingin dihapus
kolom_hapus = ['FolderPath',  'BeginTime', 'HasLabel', 'LabelID']

# Hapus kolom
loc_spring2 = loc_spring2.drop(columns=kolom_hapus)
loc_spring2.head()

# Simpan shapefile yang telah diubah ke dalam file baru
loc_spring2.to_file('loc_spring2.shp')
In [7]:
wilpen = gpd.read_file("wilpenfinal2.shp")

wilpen.head(2)
Out[7]:
KODE_UNSUR NAMA_UNSUR DESA KECAMATAN KABUPATEN PROVINSI PELAKSANA UPDATED luas geometry
0 40204 Desa Bantarmara Village KECAMATAN CISARUA KABUPATEN SUMEDANG NaN PT. DIGITAL IMAGING GEOSPATIAL 20161003 152 POLYGON Z ((830126.700 9243263.979 0.000, 8301...
1 40204 Desa Cisalak Village KECAMATAN CISARUA KABUPATEN SUMEDANG NaN PT. DIGITAL IMAGING GEOSPATIAL 20161003 231 POLYGON Z ((825787.857 9244365.502 0.000, 8257...
In [8]:
litologi = gpd.read_file("geolfinal.shp")

litologi.head(2)
Out[8]:
FID_newaja Id FID_geol Id_1 namageol geol litologi waktu Shape_Leng Shape_Area area_m Jenis area_ha geometry
0 0 0 1 0 Hasil gunung api muda, tak teruraikan Hasil gunung api muda Lava andesit Kuarter 0.000000 0.000000e+00 1.162219e+07 Qyu 1162.220 MULTIPOLYGON Z (((828365.769 9242548.483 0.000...
1 0 0 2 0 Formasi Kaliwangu Formasi Kaliwangu Konglomerat dan batu pasir Tersier 9993.148148 1.714135e+06 1.559182e+06 Pk 155.918 MULTIPOLYGON Z (((828281.013 9243420.180 0.000...
In [9]:
import geopandas as gpd

loc_spring2 = loc_spring2.to_crs('EPSG:4326')  # Proyeksi GCS_WGS_1984
wilpen = wilpen.to_crs('EPSG:4326')
litologi = litologi.to_crs('EPSG:4326')
In [10]:
#ambil nama desanya aja, gausah pakai kata 'village'
wilpen['DESA'] = wilpen['DESA'].str.split().str[0]

Visualisasi SHP¶

In [11]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots()

# Menampilkan shapefile pertama
wilpen.plot(ax=ax, color='blue', label='Kecamatan Cisarua')

# Menampilkan shapefile kedua
loc_spring2.plot(ax=ax, color='red', label='Lokasi Mata Air')

# legenda
ax.legend()

plt.show()

Peta Interaktif 1¶

In [12]:
#menampilkan di google map

import folium

# Membuat peta folium
m = folium.Map(location=[wilpen.geometry.centroid.y.mean(), wilpen.geometry.centroid.x.mean()], zoom_start=14)

# Menambahkan shapefile wilpen ke peta
folium.GeoJson(wilpen, name='Kecamatan Cisarua').add_to(m)

# Menambahkan shapefile wilpen ke peta
folium.GeoJson(litologi, name='litologi').add_to(m)

# Menambahkan shapefile loc_spring ke peta
folium.GeoJson(loc_spring2, name='Lokasi Mata Air').add_to(m)

# legenda
layer_control = folium.LayerControl()
layer_control.add_to(m)

# Menampilkan peta
m
Out[12]:
Make this Notebook Trusted to load map: File -> Trust Notebook

Peta Interaktif 2¶

In [13]:
import folium

# Membuat objek peta dengan lokasi dan zoom awal
m = folium.Map(location=[wilpen.centroid.y.mean(), wilpen.centroid.x.mean()], zoom_start=14.17)

# Cbikin layer
wilpen_layer = folium.FeatureGroup(name='Administrasi Kecamatan Cisarua')

# Menambahkan layer polygon wilpen ke dalam peta
folium.GeoJson(wilpen, name='Administrasi Kecamatan Cisarua',
               style_function=lambda x: {'fillColor': 'white', 'fillOpacity': 0.4,
                                         'color': 'black', 'weight': 1},
               highlight_function=lambda x: {'fillColor': 'yellow', 'fillOpacity': 0.6},
               tooltip=folium.GeoJsonTooltip(fields=['DESA', 'luas'],
                                              aliases=['Desa/Kelurahan', 'Luas (m)'])).add_to(wilpen_layer)

# masukin
wilpen_layer.add_to(m)

# bikin layer
litologi_layer = folium.FeatureGroup(name='Litologi')


# Menambahkan layer polygon litologi ke dalam peta
folium.GeoJson(litologi, name='Litologi',
               style_function=lambda x: {'fillColor': 'pink', 'fillOpacity': 0.4,
                                         'color': 'black', 'weight': 1},
               highlight_function=lambda x: {'fillColor': 'yellow', 'fillOpacity': 0.6},
               tooltip=folium.GeoJsonTooltip(fields=['litologi', 'geol'],
                                              aliases=['Litologi', 'Geologi'])).add_to(litologi_layer)

# masukin
litologi_layer.add_to(m)


# bikin layer
loc_spring2_layer = folium.FeatureGroup(name='Lokasi Mata Air')

# Menambahkan layer point loc_spring2 ke dalam layer grup
for idx, row in loc_spring2.iterrows():
    folium.CircleMarker(location=[row['geometry'].y, row['geometry'].x],
                        radius=row['debitt2'],
                        color='green',
                        fill=False,
                        fill_color='green',
                        fill_opacity=1,
                        tooltip=row['Name'],
                        popup=f"PH: {row['PH']}, Debit: {row['debitt2']}").add_to(loc_spring2_layer)


#masukin
loc_spring2_layer.add_to(m)

# Menambahkan legenda untuk masing-masing layer
folium.LayerControl().add_to(m)

# Menambahkan judul peta
title_html = '''
             <h3 align="center" style="font-size:30px"><b>Lokasi Mata Air di Kabupaten Sumedang</b></h3>
             '''
m.get_root().html.add_child(folium.Element(title_html))

# Menampilkan peta
m
Out[13]:
Make this Notebook Trusted to load map: File -> Trust Notebook
In [14]:
!pip install geopandas plotly mapboxgl
Requirement already satisfied: geopandas in c:\users\user\anaconda3\lib\site-packages (0.13.2)
Requirement already satisfied: plotly in c:\users\user\anaconda3\lib\site-packages (5.6.0)
Requirement already satisfied: mapboxgl in c:\users\user\anaconda3\lib\site-packages (0.10.2)
Requirement already satisfied: fiona>=1.8.19 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (1.9.4.post1)
Requirement already satisfied: shapely>=1.7.1 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (2.0.1)
Requirement already satisfied: packaging in c:\users\user\anaconda3\lib\site-packages (from geopandas) (21.3)
Requirement already satisfied: pyproj>=3.0.1 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (3.6.0)
Requirement already satisfied: pandas>=1.1.0 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (1.4.2)
Requirement already satisfied: tenacity>=6.2.0 in c:\users\user\anaconda3\lib\site-packages (from plotly) (8.0.1)
Requirement already satisfied: six in c:\users\user\anaconda3\lib\site-packages (from plotly) (1.16.0)
Requirement already satisfied: geojson in c:\users\user\anaconda3\lib\site-packages (from mapboxgl) (3.0.1)
Requirement already satisfied: matplotlib in c:\users\user\anaconda3\lib\site-packages (from mapboxgl) (3.5.1)
Requirement already satisfied: jinja2 in c:\users\user\anaconda3\lib\site-packages (from mapboxgl) (2.11.3)
Requirement already satisfied: colour in c:\users\user\anaconda3\lib\site-packages (from mapboxgl) (0.1.5)
Requirement already satisfied: chroma-py in c:\users\user\anaconda3\lib\site-packages (from mapboxgl) (0.1.0.dev1)
Requirement already satisfied: cligj>=0.5 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: click~=8.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (8.0.4)
Requirement already satisfied: attrs>=19.2.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (21.4.0)
Requirement already satisfied: importlib-metadata in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (4.11.3)
Requirement already satisfied: click-plugins>=1.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: certifi in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (2023.5.7)
Requirement already satisfied: colorama in c:\users\user\anaconda3\lib\site-packages (from click~=8.0->fiona>=1.8.19->geopandas) (0.4.4)
Requirement already satisfied: numpy>=1.18.5 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (1.22.4)
Requirement already satisfied: pytz>=2020.1 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2021.3)
Requirement already satisfied: python-dateutil>=2.8.1 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2.8.2)
Requirement already satisfied: zipp>=0.5 in c:\users\user\anaconda3\lib\site-packages (from importlib-metadata->fiona>=1.8.19->geopandas) (3.7.0)
Requirement already satisfied: MarkupSafe>=0.23 in c:\users\user\anaconda3\lib\site-packages (from jinja2->mapboxgl) (2.0.1)
Requirement already satisfied: fonttools>=4.22.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->mapboxgl) (4.25.0)
Requirement already satisfied: cycler>=0.10 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->mapboxgl) (0.11.0)
Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->mapboxgl) (1.3.2)
Requirement already satisfied: pyparsing>=2.2.1 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->mapboxgl) (3.0.4)
Requirement already satisfied: pillow>=6.2.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->mapboxgl) (9.0.1)

Peta Interaktif 3¶

In [15]:
#menampilkan di mapbox

import geopandas as gpd
import shapely.geometry as sg
import plotly.graph_objects as go
import plotly.express as px
import plotly.io as pio
import mapboxgl

# Membuat peta menggunakan plotly express
fig = px.choropleth_mapbox(wilpen, geojson=wilpen.geometry.__geo_interface__,
                           locations=wilpen.index, color="DESA",
                           color_discrete_sequence=px.colors.cyclical.IceFire, opacity=0.4,
                           center={"lat": wilpen.centroid.y.mean(), "lon": wilpen.centroid.x.mean()},
                           hover_name="DESA",
                           mapbox_style="carto-positron", zoom=13)

# Menambahkan trace untuk shapefile point loc_spring2
scatter_trace = px.scatter_mapbox(
    loc_spring2,
    lat=loc_spring2.geometry.y,
    lon=loc_spring2.geometry.x,
    color=loc_spring2['PH'],
    size=loc_spring2['debitt2'],
    opacity=1,
    text=loc_spring2['Name'],
    hover_name="Name"
).data[0]

# masukin trace
fig = fig.add_trace(scatter_trace)

# Update colorbar properties
fig.update_layout(
    coloraxis=dict(
        colorbar=dict(
            title='Nilai PH',
            x=1.025,
            y=0.3,
            len=0.5,
            tickfont=dict(size=10),
            tickformat=".1f",
            title_font=dict(size=12)
        )
    )
)

# Set the legend title and itemsizing
fig.update_layout(
    legend={"title": "Desa/Kelurahan", "itemsizing": "constant"},
    showlegend=True)

# Set the map title
fig.update_layout(
    title='Lokasi Mata Air di Kecamatan Cisarua',
    title_font=dict(size=24),
    coloraxis_colorbar=dict(title='Nilai PH', len=0.5, tickformat='.1f'),
    coloraxis=dict(colorbar=dict(y=0.3)))

fig.show()

print('*Makin besar ukuran point lokasi Mata Air, maka makin besar pula debitnya')
*Makin besar ukuran point lokasi Mata Air, maka makin besar pula debitnya

Peta Interaktif 3¶

In [16]:
fig2 = px.choropleth_mapbox(litologi, geojson=litologi.geometry.__geo_interface__,
                           locations=litologi.index, color="litologi",
                           color_discrete_sequence=px.colors.cyclical.Twilight, opacity=0.4,
                           center={"lat": litologi.centroid.y.mean(), "lon": litologi.centroid.x.mean()},
                           hover_name="litologi",
                           mapbox_style="carto-positron", zoom=12.5)

# Menambahkan trace untuk shapefile point loc_spring2
scatter_trace = px.scatter_mapbox(
    loc_spring2,
    lat=loc_spring2.geometry.y,
    lon=loc_spring2.geometry.x,
    color=loc_spring2['PH'],
    size=loc_spring2['debitt2'],
    opacity=1,
    text=loc_spring2['Name'],
    hover_name="Name"
).data[0]

fig2.add_trace(scatter_trace)

# Set the map title
fig2.update_layout(
    title='Lokasi Mata Air di Kecamatan Cisarua Berdasarkan Litologinya',
    title_font=dict(size=24),
    coloraxis_colorbar=dict(title='Nilai PH', len=0.5, tickformat='.1f'),
    coloraxis=dict(colorbar=dict(y=0.3))
)


# Set the legend
fig2.update_layout(
    legend=dict(
        title='Litologi',
        itemsizing='constant'
    )
)

fig2.show()

print('*Makin besar ukuran point lokasi Mata Air, maka makin besar pula debitnya')
*Makin besar ukuran point lokasi Mata Air, maka makin besar pula debitnya
In [17]:
# cek warna
import plotly.express.colors as px_colors

available_colors_sequential = [name for name in dir(px_colors.sequential) if not name.startswith('__')]
available_colors_diverging = [name for name in dir(px_colors.diverging) if not name.startswith('__')]

print("Sequential Colors:")
print(available_colors_sequential)

print("Diverging Colors:")
print(available_colors_diverging)
Sequential Colors:
['Aggrnyl', 'Aggrnyl_r', 'Agsunset', 'Agsunset_r', 'Blackbody', 'Blackbody_r', 'Bluered', 'Bluered_r', 'Blues', 'Blues_r', 'Blugrn', 'Blugrn_r', 'Bluyl', 'Bluyl_r', 'Brwnyl', 'Brwnyl_r', 'BuGn', 'BuGn_r', 'BuPu', 'BuPu_r', 'Burg', 'Burg_r', 'Burgyl', 'Burgyl_r', 'Cividis', 'Cividis_r', 'Darkmint', 'Darkmint_r', 'Electric', 'Electric_r', 'Emrld', 'Emrld_r', 'GnBu', 'GnBu_r', 'Greens', 'Greens_r', 'Greys', 'Greys_r', 'Hot', 'Hot_r', 'Inferno', 'Inferno_r', 'Jet', 'Jet_r', 'Magenta', 'Magenta_r', 'Magma', 'Magma_r', 'Mint', 'Mint_r', 'OrRd', 'OrRd_r', 'Oranges', 'Oranges_r', 'Oryel', 'Oryel_r', 'Peach', 'Peach_r', 'Pinkyl', 'Pinkyl_r', 'Plasma', 'Plasma_r', 'Plotly3', 'Plotly3_r', 'PuBu', 'PuBuGn', 'PuBuGn_r', 'PuBu_r', 'PuRd', 'PuRd_r', 'Purp', 'Purp_r', 'Purples', 'Purples_r', 'Purpor', 'Purpor_r', 'Rainbow', 'Rainbow_r', 'RdBu', 'RdBu_r', 'RdPu', 'RdPu_r', 'Redor', 'Redor_r', 'Reds', 'Reds_r', 'Sunset', 'Sunset_r', 'Sunsetdark', 'Sunsetdark_r', 'Teal', 'Teal_r', 'Tealgrn', 'Tealgrn_r', 'Turbo', 'Turbo_r', 'Viridis', 'Viridis_r', 'YlGn', 'YlGnBu', 'YlGnBu_r', 'YlGn_r', 'YlOrBr', 'YlOrBr_r', 'YlOrRd', 'YlOrRd_r', '_cols', '_contents', '_k', '_swatches', '_swatches_continuous', 'algae', 'algae_r', 'amp', 'amp_r', 'deep', 'deep_r', 'dense', 'dense_r', 'gray', 'gray_r', 'haline', 'haline_r', 'ice', 'ice_r', 'matter', 'matter_r', 'solar', 'solar_r', 'speed', 'speed_r', 'swatches', 'swatches_continuous', 'tempo', 'tempo_r', 'thermal', 'thermal_r', 'turbid', 'turbid_r']
Diverging Colors:
['Armyrose', 'Armyrose_r', 'BrBG', 'BrBG_r', 'Earth', 'Earth_r', 'Fall', 'Fall_r', 'Geyser', 'Geyser_r', 'PRGn', 'PRGn_r', 'PiYG', 'PiYG_r', 'Picnic', 'Picnic_r', 'Portland', 'Portland_r', 'PuOr', 'PuOr_r', 'RdBu', 'RdBu_r', 'RdGy', 'RdGy_r', 'RdYlBu', 'RdYlBu_r', 'RdYlGn', 'RdYlGn_r', 'Spectral', 'Spectral_r', 'Tealrose', 'Tealrose_r', 'Temps', 'Temps_r', 'Tropic', 'Tropic_r', '_cols', '_contents', '_k', '_swatches', '_swatches_continuous', 'balance', 'balance_r', 'curl', 'curl_r', 'delta', 'delta_r', 'oxy', 'oxy_r', 'swatches', 'swatches_continuous']
In [18]:
#palet warna
import plotly.express.colors as px_colors

cyclical_palettes = [name for name in dir(px_colors.cyclical) if not name.startswith('__')]

print(cyclical_palettes)
['Edge', 'Edge_r', 'HSV', 'HSV_r', 'IceFire', 'IceFire_r', 'Phase', 'Phase_r', 'Twilight', 'Twilight_r', '_cols', '_contents', '_k', '_swatches', '_swatches_continuous', '_swatches_cyclical', 'mrybm', 'mrybm_r', 'mygbm', 'mygbm_r', 'swatches', 'swatches_continuous', 'swatches_cyclical']

Visualisasi Interpolasi¶

In [19]:
!pip install geopandas pykrige numpy matplotlib
Requirement already satisfied: geopandas in c:\users\user\anaconda3\lib\site-packages (0.13.2)
Requirement already satisfied: pykrige in c:\users\user\anaconda3\lib\site-packages (1.7.0)
Requirement already satisfied: numpy in c:\users\user\anaconda3\lib\site-packages (1.22.4)
Requirement already satisfied: matplotlib in c:\users\user\anaconda3\lib\site-packages (3.5.1)
Requirement already satisfied: packaging in c:\users\user\anaconda3\lib\site-packages (from geopandas) (21.3)
Requirement already satisfied: pandas>=1.1.0 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (1.4.2)
Requirement already satisfied: shapely>=1.7.1 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (2.0.1)
Requirement already satisfied: pyproj>=3.0.1 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (3.6.0)
Requirement already satisfied: fiona>=1.8.19 in c:\users\user\anaconda3\lib\site-packages (from geopandas) (1.9.4.post1)
Requirement already satisfied: scipy<2,>=1.1.0 in c:\users\user\anaconda3\lib\site-packages (from pykrige) (1.7.3)
Requirement already satisfied: pyparsing>=2.2.1 in c:\users\user\anaconda3\lib\site-packages (from matplotlib) (3.0.4)
Requirement already satisfied: cycler>=0.10 in c:\users\user\anaconda3\lib\site-packages (from matplotlib) (0.11.0)
Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\user\anaconda3\lib\site-packages (from matplotlib) (1.3.2)
Requirement already satisfied: pillow>=6.2.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib) (9.0.1)
Requirement already satisfied: python-dateutil>=2.7 in c:\users\user\anaconda3\lib\site-packages (from matplotlib) (2.8.2)
Requirement already satisfied: fonttools>=4.22.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib) (4.25.0)
Requirement already satisfied: importlib-metadata in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (4.11.3)
Requirement already satisfied: click-plugins>=1.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.1.1)
Requirement already satisfied: certifi in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (2023.5.7)
Requirement already satisfied: cligj>=0.5 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (0.7.2)
Requirement already satisfied: attrs>=19.2.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (21.4.0)
Requirement already satisfied: six in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (1.16.0)
Requirement already satisfied: click~=8.0 in c:\users\user\anaconda3\lib\site-packages (from fiona>=1.8.19->geopandas) (8.0.4)
Requirement already satisfied: colorama in c:\users\user\anaconda3\lib\site-packages (from click~=8.0->fiona>=1.8.19->geopandas) (0.4.4)
Requirement already satisfied: pytz>=2020.1 in c:\users\user\anaconda3\lib\site-packages (from pandas>=1.1.0->geopandas) (2021.3)
Requirement already satisfied: zipp>=0.5 in c:\users\user\anaconda3\lib\site-packages (from importlib-metadata->fiona>=1.8.19->geopandas) (3.7.0)
In [20]:
import geopandas as gpd
from scipy.interpolate import Rbf
import numpy as np
import matplotlib.pyplot as plt
In [21]:
import geopandas as gpd
from scipy.interpolate import Rbf
import numpy as np
import matplotlib.pyplot as plt

# Extract koordinat shp 'wilpen'
x = []
y = []
for polygon in wilpen.geometry:
    if polygon.geom_type == 'Polygon':
        coords = polygon.exterior.coords
        x += [coord[0] for coord in coords]
        y += [coord[1] for coord in coords]
    elif polygon.geom_type == 'MultiPolygon':
        for subpolygon in polygon:
            coords = subpolygon.exterior.coords
            x += [coord[0] for coord in coords]
            y += [coord[1] for coord in coords]

# Memotong data titik observasi dengan polygon batas administrasi
obs_points_cropped = gpd.overlay(loc_spring2, wilpen, how='intersection')

# Ambil koordinat titik observasi yang dipotong
x_obs = obs_points_cropped['lat'].values
y_obs = obs_points_cropped['long'].values
values = obs_points_cropped['PH'].values

# Menentukan grid untuk interpolasi
grid_x = np.linspace(min(x), max(x), 100)
grid_y = np.linspace(min(y), max(y), 100)
grid_xx, grid_yy = np.meshgrid(grid_x, grid_y)
In [22]:
obs_points_cropped.columns
Out[22]:
Index(['OID_', 'Name', 'EndTime', 'kelas_pema', 'luas_jangk', 'Nama',
       'TDS__ppm_2', 'debitt2', 'lat', 'long', 'Kode', 'Mata air', 'PH',
       'Salinitas', 'Suhu', 'TDS', 'Debit', 'DHL', 'warna', 'Rasa',
       'Jarak patahan', 'kelompok', 'KODE_UNSUR', 'NAMA_UNSUR', 'DESA',
       'KECAMATAN', 'KABUPATEN', 'PROVINSI', 'PELAKSANA', 'UPDATED', 'luas',
       'geometry'],
      dtype='object')
In [23]:
# Melakukan interpolasi dengan IDW
rbf = Rbf(x_obs, y_obs, values, function='linear')
interp_values = rbf(grid_xx, grid_yy)
In [24]:
interp_values
Out[24]:
array([[ 9.16507837,  9.05443297,  8.94737069, ..., 10.72502986,
        10.83152736, 10.93963369],
       [ 9.10938659,  8.99741811,  8.88904567, ..., 10.67515315,
        10.78208318, 10.89063202],
       [ 9.05438672,  8.94104327,  8.8313002 , ..., 10.62561109,
        10.73297816, 10.84197413],
       ...,
       [ 9.4572802 ,  9.34707808,  9.23843677, ...,  9.18753308,
         9.32312406,  9.46034547],
       [ 9.50105329,  9.39136294,  9.28324302, ...,  9.23151407,
         9.36587144,  9.50191993],
       [ 9.54540992,  9.43623481,  9.32863917, ...,  9.27640647,
         9.40953893,  9.54441876]])
In [25]:
!pip install rasterio
Requirement already satisfied: rasterio in c:\users\user\anaconda3\lib\site-packages (1.3.8)
Requirement already satisfied: numpy>=1.18 in c:\users\user\anaconda3\lib\site-packages (from rasterio) (1.22.4)
Requirement already satisfied: cligj>=0.5 in c:\users\user\anaconda3\lib\site-packages (from rasterio) (0.7.2)
Requirement already satisfied: setuptools in c:\users\user\anaconda3\lib\site-packages (from rasterio) (61.2.0)
Requirement already satisfied: affine in c:\users\user\anaconda3\lib\site-packages (from rasterio) (2.4.0)
Requirement already satisfied: certifi in c:\users\user\anaconda3\lib\site-packages (from rasterio) (2023.5.7)
Requirement already satisfied: snuggs>=1.4.1 in c:\users\user\anaconda3\lib\site-packages (from rasterio) (1.4.7)
Requirement already satisfied: click-plugins in c:\users\user\anaconda3\lib\site-packages (from rasterio) (1.1.1)
Requirement already satisfied: click>=4.0 in c:\users\user\anaconda3\lib\site-packages (from rasterio) (8.0.4)
Requirement already satisfied: attrs in c:\users\user\anaconda3\lib\site-packages (from rasterio) (21.4.0)
Requirement already satisfied: colorama in c:\users\user\anaconda3\lib\site-packages (from click>=4.0->rasterio) (0.4.4)
Requirement already satisfied: pyparsing>=2.1.6 in c:\users\user\anaconda3\lib\site-packages (from snuggs>=1.4.1->rasterio) (3.0.4)
In [26]:
import geopandas as gpd
import numpy as np
import plotly.graph_objects as go
import rasterio
from rasterio.transform import from_origin
In [27]:
import folium
from folium import plugins
from io import BytesIO
import base64
In [ ]:
 
In [28]:
import geopandas as gpd
import matplotlib.pyplot as plt

# ambil raster
interpolated_raster = interp_values

# Create the plot
fig, ax = plt.subplots(figsize=(20, 20))
img = ax.imshow(interpolated_raster, extent=[grid_x.min(), grid_x.max(), grid_y.min(), grid_y.max()], origin='lower', cmap='coolwarm')
wilpen.plot(ax=ax, color='none', edgecolor='red', linewidth=0.9, label='Batas Administrasi Desa')
loc_spring2.plot(ax=ax, color='red', markersize=5, label='Lokasi Mata Air')

# Add a colorbar
cbar = plt.colorbar(img, ax=ax, label='Nilai PH', shrink=0.4)

# Add annotations for the shapefile items
for idx, row in wilpen.iterrows():
    label = row['DESA']
    x, y = row.geometry.centroid.x, row.geometry.centroid.y
    ax.annotate(label, (x, y), color='red', fontsize=17, ha='center', va='center')

loc_spring2_annotations = []
for idx, row in loc_spring2.iterrows():
    label = row['Nama']
    x, y = row.geometry.x, row.geometry.y
    ann = ax.annotate(label, xy=(x, y), xytext=(5, -5), textcoords='offset points',
                      ha='left', va='top', color='white', fontsize=10,
                      bbox=dict(facecolor='none', edgecolor='none', boxstyle='round,pad=0.3'))
    loc_spring2_annotations.append(ann)

# Adjust the positions of loc_spring2 annotations to avoid overlapping
plt.draw()
for i, ann in enumerate(loc_spring2_annotations):
    for j in range(i + 1, len(loc_spring2_annotations)):
        ann_j = loc_spring2_annotations[j]
        if ann.get_bbox_patch().get_extents().overlaps(ann_j.get_bbox_patch().get_extents()):
            x, _ = ann.xy
            x_j, _ = ann_j.xy
            offset = abs(x_j - x) + 2
            ann_j.xy = (x_j + offset, ann_j.xy[1])

ax.set_title('Prediksi PH Air Tanah Melalui Mata Air\ndi Kecamatan Cisarua, Sumedang\n(menggunakan IDW)', fontsize = 25, pad=15)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

plt.show()
In [29]:
# Extract koordinat dr shp litologi
x = []
y = []

for geometry in litologi.geometry:
    if geometry.geom_type == 'Polygon':
        coords = geometry.exterior.coords
        x += [coord[0] for coord in coords]
        y += [coord[1] for coord in coords]
    elif geometry.geom_type == 'MultiPolygon':
        for polygon in geometry.geoms:
            coords = polygon.exterior.coords
            x += [coord[0] for coord in coords]
            y += [coord[1] for coord in coords]

# Memotong data titik observasi dengan polygon batas administrasi
obs_points_cropped = gpd.overlay(loc_spring2, litologi, how='intersection')
In [30]:
obs_points_cropped.columns
Out[30]:
Index(['OID_', 'Name', 'EndTime', 'kelas_pema', 'luas_jangk', 'Nama',
       'TDS__ppm_2', 'debitt2', 'lat', 'long', 'Kode', 'Mata air', 'PH',
       'Salinitas', 'Suhu', 'TDS', 'Debit', 'DHL', 'warna', 'Rasa',
       'Jarak patahan', 'kelompok', 'FID_newaja', 'Id', 'FID_geol', 'Id_1',
       'namageol', 'geol', 'litologi', 'waktu', 'Shape_Leng', 'Shape_Area',
       'area_m', 'Jenis', 'area_ha', 'geometry'],
      dtype='object')
In [31]:
# Ambil koordinat titik observasi yang dipotong
x_obs = obs_points_cropped['lat'].values
y_obs = obs_points_cropped['long'].values
values = obs_points_cropped['PH'].values

# Menentukan grid untuk interpolasi
grid_x = np.linspace(min(x), max(x), 100)
grid_y = np.linspace(min(y), max(y), 100)
grid_xx, grid_yy = np.meshgrid(grid_x, grid_y)
In [32]:
!pip install adjustText
Requirement already satisfied: adjustText in c:\users\user\anaconda3\lib\site-packages (0.8)
Requirement already satisfied: numpy in c:\users\user\anaconda3\lib\site-packages (from adjustText) (1.22.4)
Requirement already satisfied: matplotlib in c:\users\user\anaconda3\lib\site-packages (from adjustText) (3.5.1)
Requirement already satisfied: packaging>=20.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (21.3)
Requirement already satisfied: kiwisolver>=1.0.1 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (1.3.2)
Requirement already satisfied: cycler>=0.10 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (0.11.0)
Requirement already satisfied: python-dateutil>=2.7 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (2.8.2)
Requirement already satisfied: fonttools>=4.22.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (4.25.0)
Requirement already satisfied: pyparsing>=2.2.1 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (3.0.4)
Requirement already satisfied: pillow>=6.2.0 in c:\users\user\anaconda3\lib\site-packages (from matplotlib->adjustText) (9.0.1)
Requirement already satisfied: six>=1.5 in c:\users\user\anaconda3\lib\site-packages (from python-dateutil>=2.7->matplotlib->adjustText) (1.16.0)
In [33]:
import matplotlib.pyplot as plt
from adjustText import adjust_text

#ambil data raster
interpolated_raster = interp_values

# Create the plot
fig, ax = plt.subplots(figsize=(20, 20))
img = ax.imshow(interpolated_raster, extent=[grid_x.min(), grid_x.max(), grid_y.min(), grid_y.max()], origin='lower', cmap='coolwarm')
litologi.plot(ax=ax, color='none', edgecolor='red', linewidth=0.9, label='Litologi')
loc_spring2.plot(ax=ax, color='red', markersize=5, label='Lokasi Mata Air')

cbar = plt.colorbar(img, ax=ax, label='Nilai PH', shrink=0.4)

# Add annotations for the shapefile items
texts = []
for idx, row in litologi.iterrows():
    label = row['litologi']
    x, y = row.geometry.centroid.x, row.geometry.centroid.y
    text = ax.annotate(label, (x, y), color='black', fontsize=17, ha='center', va='center')
    texts.append(text)

loc_spring2_annotations = []
for idx, row in loc_spring2.iterrows():
    label = row['Nama']
    x, y = row.geometry.x, row.geometry.y
    ann = ax.annotate(label, xy=(x, y), xytext=(5, -5), textcoords='offset points',
                      ha='left', va='top', color='white', fontsize=10,
                      bbox=dict(facecolor='none', edgecolor='none', boxstyle='round,pad=0.3'))
    loc_spring2_annotations.append(ann)

# Adjust the positions of loc_spring2 annotations to avoid overlapping
plt.draw()
for i, ann in enumerate(loc_spring2_annotations):
    for j in range(i + 1, len(loc_spring2_annotations)):
        ann_j = loc_spring2_annotations[j]
        if ann.get_bbox_patch().get_extents().overlaps(ann_j.get_bbox_patch().get_extents()):
            x, _ = ann.xy
            x_j, _ = ann_j.xy
            offset = abs(x_j - x) + 2
            ann_j.xy = (x_j + offset, ann_j.xy[1])

adjust_text(texts, arrowprops=dict(arrowstyle='-', color='black'), autoalign='xy', force_points=0.3, force_text=0.8)

ax.set_title('Prediksi PH Air Tanah Melalui Mata Air\ndi Kecamatan Cisarua, Sumedang\nBerdasarkan Kondisi Litologi (menggunakan IDW)', fontsize=25, pad=15)
ax.set_xlabel('Longitude')
ax.set_ylabel('Latitude')

plt.show()
In [35]:
# Simpan hasil interpolasi ke dalam file raster menggunakan rasterio
with rasterio.open('interp_raster.tif', 'w', driver='GTiff', height=interp_values.shape[0], width=interp_values.shape[1], count=1, dtype=interp_values.dtype) as dst:
    dst.write(interp_values, 1)
In [36]:
#nyoba pake metode lain
import geopandas as gpd
from pykrige.ok import OrdinaryKriging
from scipy.interpolate import Rbf

# Reproject shapefiles to UTM 48S
loc_spring2 = loc_spring2.to_crs('EPSG:32748')
wilpen = wilpen.to_crs('EPSG:32748')

# Extract coordinates and pH values from loc_spring2
x = loc_spring2.geometry.x
y = loc_spring2.geometry.y
z = loc_spring2['PH']

# Extract mask polygon from wilpen
mask_polygon = wilpen.geometry.unary_union
In [37]:
# Create grid points within the mask polygon
mask_bounds = mask_polygon.bounds
grid_x = np.linspace(mask_bounds[0], mask_bounds[2], num=100)
grid_y = np.linspace(mask_bounds[1], mask_bounds[3], num=100)
grid_xx, grid_yy = np.meshgrid(grid_x, grid_y)
grid_points = np.column_stack((grid_xx.flatten(), grid_yy.flatten()))

# Perform IDW interpolation
idw_interpolation = Rbf(x, y, z, function='linear')
idw_values = idw_interpolation(grid_points[:, 0], grid_points[:, 1])
In [38]:
# Create GeoDataFrame from the interpolated values
interpolated_data = gpd.GeoDataFrame(geometry=gpd.points_from_xy(grid_points[:, 0], grid_points[:, 1]))
interpolated_data['IDW_PH'] = idw_values

# Apply the mask polygon
interpolated_data = interpolated_data[interpolated_data.within(mask_polygon)]
In [39]:
idw_values
Out[39]:
array([9.09445383, 8.98441761, 8.87819642, ..., 9.27671163, 9.40995203,
       9.54495274])
In [40]:
# Set the output file path
output_file = 'interpolated_idw.tif'

# Set the CRS and transform parameters for the raster
transform = from_origin(grid_x[0], grid_y[0], grid_x[1] - grid_x[0], grid_y[1] - grid_y[0])

# Create the raster dataset
interpolated_dataset = rasterio.open(
    output_file,
    'w',
    driver='GTiff',
    height=grid_y.shape[0],
    width=grid_x.shape[0],
    count=1,
    dtype=idw_values.dtype,
    transform=transform
) 

# Write the interpolated values to the raster dataset
interpolated_dataset.write(idw_values.reshape((1, grid_y.shape[0], grid_x.shape[0])))

# Close the raster dataset
interpolated_dataset.close()
In [45]:
import rasterio
import plotly.graph_objects as go

# Read the interpolated raster data
interpolated_dataset = rasterio.open(output_file)
interpolated_data = interpolated_dataset.read(1)

# Get the geospatial information from the raster dataset
crs = interpolated_dataset.crs
transform = interpolated_dataset.transform

# Create the Mapbox figure
fig = go.Figure()

# Compute the grid coordinates
y_coords, x_coords = rasterio.transform.xy(transform, range(interpolated_data.shape[0]), range(interpolated_data.shape[1]))

# Add the raster layer to the figure
fig.add_trace(go.Heatmap(z=interpolated_data, colorscale='Viridis', zmin=interpolated_data.min(), zmax=interpolated_data.max(), zauto=False,
                        showscale=True, opacity=0.7, hoverinfo='skip',
                        x=x_coords, y=y_coords))

# Set the layout of the figure
fig.update_layout(mapbox_style="open-street-map", mapbox_zoom=12,
                  mapbox_center={"lat": sum(y_coords) / len(y_coords),
                                 "lon": sum(x_coords) / len(x_coords)})

fig.update_layout(title_text="Interpolasi PH Air Tanah Berdasarkan Mata Air")

# Show the figure
fig.show()
In [42]:
# Simpan hasil interpolasi ke dalam file raster menggunakan rasterio
#with rasterio.open('interp_data.tif', 'w', driver='GTiff', height=interpolated_data.shape[0], width=interpolated_data.shape[1], count=1, dtype=interpolated_data.dtype) as dst:
#    dst.write(interpolated_data, 1)
In [43]:
# Interpolasi menggunakan metode IDW
#idw_interpolation = Rbf(x, y, z, function='linear')
#interp_values = rbf(grid_xx, grid_yy)
In [44]:
# Interpolation using Kriging
#OK = OrdinaryKriging(x, y, z, variogram_model='linear')
#kriging_interpolation, kriging_variance = OK.execute('grid', x, y)